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Table 1: Constraints to the parameters of expanded N^j . Blanks mean no constraints. 
X0...X3 are not constraints, they are used later to calculate other parameters. 
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Figure 1: Sound speed as a function of internal 
energy e. 
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Figure 2: Shear viscosity as a function of re- 
laxation time r. e = 0.4,0.5,0.9,1.0, star. 
TLBE; circle: ETLBE; line: Theory. 
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Figure 3: Thermal conductivity as a function 
of relaxation time r. e = 0.4,0.5,0.55,0.75, 
star. TLBE; circle: ETLBE; line: Theory. 



Figure 4: Ma effects on the decaying rate 
of total kinetic energy in shear wave flow. 
V2Ma = 0.0158,0.63,0.95, dots: TLBE; 
lines: ETLBE. 
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Abstract 

We present a new thermal lattice BGK model in D-dimensional space for the numerical 
calculation of fluid dynamics. This model uses a higher order expansion of equilibrium 
distribution in Maxwellian type. In the mean time the lattice symmetry is upgraded to 
ensure the isotropy of 6th order tensor. These manipulations lead to macroscopic equations 
free from nonlinear deviations. We demonstrate the improvements by conducting classical 
Chapman-Enskog analysis and by numerical simulation of shear wave flow. The transport 
coefficients are measured numerically, too. 



1 INTRODUCTION 

Since the proposal of lattice gas automata(LGA) 
model for the simulation of fluid flow, much at- 
tention has been paid to the totally discrete ki- 
netic theory. The governing equation for this rela- 
tively new theory is called as the lattice Boltzmann 
equation(LBE). Although being a derivation of the 
LGA model, LBE became popular rapidly among 
the enthusiastic researchers of lattice hydrodynam- 
ics because it is noise-free and so flexible that the 
physics of real fluid could be fully recovered. These 
advantages may well be observed from the newly 
suggested lattice BGK modelfl], in which the BGK 
collision operator[2] was introduced to the origi- 
nal LBE models. Today, main reasons that could 
be forwarded to persuade a fluid dynamic expert 
to adopt the lattice BGK model as an alternative 
tool for numerical simulation lie in three folds: the 
simple algorithm, the easiness in dealing with com- 
plicated geometric boundaries and the high level of 
parallelism in the implementation. 

The recent topics about the lattice BGK mod- 
els concerned with the nonlinear deviation in com- 
pressible regime and models that consider thermal 
effects. The nonlinear deviation was first derived 
by Qian[3] to be an additional dissipative term in 
the r. h. s. of the macroscopic momentum equa- 



tion of the non-thermal lattice BGK model. Hence 
the macro-dynamics of such model is actually gov- 
erned by 

Navier-Stokes Equation + adpd^(pu a upu^) . (1) 

Here p, u are density and flow velocity, a is de- 
fined as — r), where r is the relaxation time 
of the BGK collision operator. When the flow en- 
ters into the compressible regime, the deviation be- 
comes significant as the Mach number is no longer 
small. The pioneer work of the thermal lattice 
BGK model was published in [4], where the hex- 
ongal lattice and 2nd order expansion of equilib- 
rium distribution are employed. The conservation 
of mass, momentum and energy in micro- world 
gives rise to macroscopic equations which look like 
those of real fluid, except for some hidden nonlin- 
ear terms. 

The current thermal lattice BGK models thus 
are not free from nonlinear deviation. The de- 
viation in momentum equation is similar to the 
term appearing in equation (1). Deviations in 
energy equation are found to be proportional to 
dadpipeUaUp), d a df}(pu a up), etc. In order 
to eliminate these deviations, higher order expan- 
sion(4th) has to be considered for the equilibrium 
distribution function. Furthermore, the lattice 



should be symmetric enough to ensure isotropy up 
to 6th order tensor so that higher order expansions 
could be included correctly. This cannot be real- 
ized on lattices consisting of regular polygon(such 
as hexongal lattice), as they ensure at most the 4th 
order tensor to be isotropic. 

In the sequel, we shall introduce a model which 
is able to eliminate the nonlinear deviations. We 
present the theoretical works concerning with the 
micro- and macro-dynamics in Section 2. In Sec- 
tion 3 the transport coefficients will be numeri- 
cally measured and the elimination of the nonlinear 
term in the momentum equation will be illustrated. 
Section 4 gives the concluding remarks. 

2 LATTICE HYDRODYNAMICS 

2.1 Lattice Geometry and Symmetry 

As in our previous model[5], the lattice we em- 
ployed is composed of several sub-lattices. The co- 
ordinate of base vectors of these sub-lattices may 
be represented by fc(±l, . . . , ±1, 0, . . . , 0) and its 
permutations. The number of non-zero compo- 
nents is recorded as p, which also stands for the 
modulus of such base vector. A; is a multiplier for 
the magnitude of modulus. Hence the vector of lat- 
tice link may be denoted as c p ki, where pk gives the 
index of sub-lattice and i counts total b p vectors 
lie in the pk sub-lattice. The tensor of n-velocity 
moments, defined as T^ a e = J2i c pkia- ■ -c p ki(, is 
crucial to the hydrodynamic derivation. The odd 
order tensors vanish naturally by the definition it- 
self and the 2nd, 4th and 6th order tensors are re- 
quired to be isotropic. It can be shown that these 
even order tensors are generally written as 
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where T and 6 are Kronecker tensors, ". . ." rep- 
resents permutations of the indices of the first 
term. The specific values of ip p j. , ip p k , A p k , Q, p k and 
& p k could be calculated by the products of gauge 
vectors(±c p ij', ic p n ± c p \j , . . .) and are listed as 
follows for two dimensional space, 
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In case k ^ 1, these values should be multiplied by 
k n . With careful inspection of these parameters, 
we may infer that anisotropic parts of 4th and 6th 
order tensors would vanish simultaneously if the 
ratio of the particle distributions on different sub- 
lattices is properly tuned. 

2.2 Lattice BGK Equation 

Lattice BGK equation describes the propagation 
and collision of particles occurring on the discrete 
spatial lattice, at discrete time step and with dis- 
crete velocity set. Taking the particle distribution 
on each lattice link as N p ki(x,t), this equation ap- 
pears to be 



N p ki(x + c pk i,t + 1) - N pk i(x,t) = 



(5) 



Here the collisions of particles are replaced by a re- 
laxation process, in which the particle distribution 
is relaxed to its equilibrium value over a time pe- 
riod r. In order to reproduce Navier-Stokes fluid 
the equilibrium distribution should be Maxwellian 
and depends only on the local conserved density 
of mass, momentum and energy, which are defined 
by the following formulas, 

^ Npki = P , ^2 Np k iCp k i a = pu a , (6) 
pki pki 

Y, N P kiC lki = />(2e + « 2 ). 

pki 

Here e is the thermal energy determined by flying 
speed of particles. When the flow velocity is very 
much smaller than this speed, the equilibrium dis- 
tribution is expanded into the Chapman-Enskog 
form with the consideration of parity invariance of 
regular lattices, 
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The parameters stated above depend only on local 
p and e and are to be written in this form, 



Xpk = p(x p ko + Xpkie + Xpk-ze 2 ) ■ 



(8) 



Xpk may represent any of A p k,M p k . . . etc. The 
specific values of these parameters could be de- 
termined by considering higer order isotropy and 
matching the form of macroscopic equations with 
those of physical fluids. We summarize constraints 
to the parameters that would lead to the N-S equa- 
tion in Table 1. 

2.3 Macroscopic Equations 

Using the constraints listed above and the sym- 
metric properties given in equation (2), we cal- 
culated the velocity moments of the equilibrium 
distribution to 4th order. The obtained results 
not only verified the definitions of mass, momen- 
tum and energy but also gave isotropic momen- 
tum and energy fluxes. Our results could be 
reduced to the ones published in [6], in which 
only the special three dimensional case was con- 
sidered. Next we employed the multi-scale tech- 
nique[5] and derived the macroscopic equations in 
a straightforward way. The resulting equations 
are exactly the compressible Navier-Stokes equa- 
tions(equations(l), (2), (3) in [4]). Here, in D di- 
mensions, we list the equation of state, 
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and give expressions for transport coefficients, 
2 ( l \ 
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k = — — -pe(T ). 

D 2 V 2 ' 

Here p is the thermodynamic pressure and p, A 
and k stand for shear viscosity, second viscosity 
and thermal conductivity respectively. 

We observed from analytical results that terms 
of nonlinear deviation no longer bother the macro- 
dynamics of the lattice BGK model. Partly 
because of this and partly because of conve- 
nience, we label our model in the sequel as 
ETLBE(Exact Thermal LBE) and the conven- 
tional model TLBE(Thermal LBE). 

3 NUMERICAL RESULTS 

We carried out all the numerical calculation on 
a two dimensional least ETLBE, in which 16 ve- 
locities are employed and indices of sub-lattices are 



11, 12, 21 and 22. To show the improvement of the 
new model, we also coded a 13-velocity TLBE in 
which equilibrium distribution is expanded to 2nd 
order and isotropy is kept only to 4th order. The 
numerical results are compared with each other for 
every case. 

3.1 Transport Coefficients 

We numerically measured the sound speed, shear 
viscosity and thermal conductivity by using the lin- 
ear perturbation theory. Figure 1 is the results of 
sound speed. The results of two models agree well 
with the theoretical prediction. The results of mea- 
surement of shear viscosity are shown in Figure 2. 
We observed deviation from theoretical value for 
both ETLBE and TLBE when r gets large enough. 
This is due to increment of Knudsen number which 
may break the hydrodynamic mode. It is interest- 
ing, however, that the way the deviation increases 
is different regarding to each model at higher en- 
ergy level. The heat conductivity is measured by 
switching off the convection effects. Both the mod- 
els give good results, see Figure 3. 

3.2 Elimination of Nonlinear Deviation 

We illustrate the elimination by simulating the 
shear wave flow 

u x = u 0) u y = u e i( " t+kx \ (11) 

on the 64x64 lattice field. The flow velocity «o 
could be changed to adjust Mach number. The 
Mach number effects on the decaying rate of the 
total kinetic energy are investigated and shown in 
Figure 4. We observed that the ETLBE behaves 
extremely well even y/2Ma reaches 0.95. The de- 
caying rates consistently stick to the value of 2vk 2 , 
where k is the wave number. For the TLBE, the 
deviation is intolerable. However this deviation 
agrees with the theoretically predicted nonlinear 
term as well, namely, the total kinetic energy de- 
cays as e 



-2(v-au' )k't 



4 CONCLUSION 

We developed a new thermal lattice BGK model 
which eliminates the reported nonlinear deviations 
effectively. This was demonstrated both analyti- 
cally and numerically. A much detailed compari- 
son between the results of lattice BGK simulation 
and the Navier-Stokes solution is very much ex- 
pected. 



